1 ARMA(p, q) Model Rewrite

Recall ARMA model and its backshift notation. We rewrite ${\displaystyle \varphi (B)(y_t-\mu )=\theta (B){\epsilon }_t}$ to $$\varphi (B)y_t=\delta +\theta (B){\epsilon }_t$$ and try to determine ${\displaystyle \delta }$. We first write $$y_t-\mu =\frac{\theta (B)}{\varphi (B)}{\epsilon }_t.$$ Factorize ${\displaystyle \varphi (z)}$ as $$\varphi (z)=(1-a_{1}z)\cdots (1-a_{p}z).$$ Here ${\displaystyle a_1^{-1},\cdots ,a_p^{-1}}$ are roots of ${\displaystyle \varphi (z)}$. Then $$y_t-\mu =\frac{\theta (B)}{\prod _{k=1}^p(1-a_{k}B)}{\epsilon }_t=\theta (B)(1-a_{1}B{)}^{-1}\cdots (1-a_{p}B{)}^{-1}{\epsilon }_t.$$
Recall expansion in (2.4), if all ${\displaystyle |a_k|<1}$, rewrite $$y_t=\mu +\sum _{j=0}^{\mathrm{\infty }}{\psi }_j{\epsilon }_{t-j}.$$ This is a causal stationary process.

Since PACF and ACF have not a clear cut off, we usually use AIC, BIC to determine ${\displaystyle p,q}$.

2 Box-Jenkins Time Series Modeling Strategy

The idea is

• Transform ${\displaystyle y_1,\cdots ,y_n}$ to ${\displaystyle x_t}$ that doesn't have any discernible trends.
• Fit an ${\displaystyle \mathrm{ARMA}(p,q)}$ model.

To implement the first idea, we usually have two ways:

• Differencing: use ${\displaystyle \mathrm{\nabla }{y}_t=y_t-y_{t-1}}$, ${\displaystyle {\mathrm{\nabla }}^2{y}_t=y_t-2y_{t-1}+y_{t-2}}$, etc.
• Seasonal Differencing: ${\displaystyle {\mathrm{\nabla }}_s{y}_t=y_t-y_{t-s}}$, ${\displaystyle {\mathrm{\nabla }}_s^2{y}_t=y_t-2y_{t-s}+y_{t-2s}}$.

3 ARIMA Models

3.1 Seasonal ARMA Models

We say ${\displaystyle \{y_t\}}$ is a seasonal ${\displaystyle \mathrm{ARMA}(P,Q)}$ with period ${\displaystyle s}$, if ${\displaystyle \mathrm{\Phi }(B^s)(y_t-\mu )=\mathrm{\Theta }(B^s){\epsilon }_t}$, where $$\begin{array}{rl}\mathrm{\Phi }(B^s)& =1-{\mathrm{\Phi }}_1{B}^s-\cdots -{\mathrm{\Phi }}_P{B}^{Ps},\\ \mathrm{\Theta }(B^s)& =1+{\mathrm{\Theta }}_1{B}^s+\cdots +{\mathrm{\Theta }}_Q{B}^{Qs}.\end{array}$$
This is a special case of an ${\displaystyle \mathrm{ARMA}(Ps,Qs)}$ model. But it has ${\displaystyle P+Q+1}$ parameters (1 for ${\displaystyle {\sigma }^2}$) while a general ${\displaystyle \mathrm{ARMA}(Ps,Qs)}$ has ${\displaystyle Ps+Qs+1}$.

The ACF and PACF of seasonal ARMA are non-zero only at seasonal lags ${\displaystyle h=0.s,2s,\cdots }$. At seasonal lags, PACF and ACF behave just like unseasonal ARMA model: ${\displaystyle \mathrm{\Phi }(B)X_t=\mathrm{\Theta }(B){\epsilon }_t}$.

3.2 Multiplicative Seasonal ARMA Models

For a dataset that might have sample autocorrelations nonnegligile at lags ${\displaystyle 0,1,11,12,13}$ (like co2 dataset), we can use this to reduce parameter.

4 SARIMA Models

5 Parameter Estimation in MA(1)

Estimating parameters of ARMA/ARIMA/SARIMA is much harder than AR models. We'll illustrate the difficulty using the example of ${\displaystyle \mathrm{MA}(1)}$. Recall from here, ${\displaystyle \mathrm{MA}(1)}$ is given by $$\begin{array}{}\text{(5.1)}& y_t=\mu +{\epsilon }_t+\theta {\epsilon }_{t-1}=\mu +\theta (B){\epsilon }_t,\theta (B)=1+\theta B,{\epsilon }_t\stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}}{\sim }\mathcal{N}(0,{\sigma }^2).\end{array}$$
The joint density of ${\displaystyle y_1,\cdots ,y_n}$ is multivariate normal with mean ${\displaystyle m=(\mu ,\cdots ,\mu {)}^{\mathrm{T}}}$ and covariance matrix ${\displaystyle \mathrm{\Sigma }}$: $$\mathrm{\Sigma }(i,j)=\{\begin{array}{rl}& {\sigma }^2(1+{\theta }^2),i=j,\\ & {\sigma }^2\theta ,|i-j|=1,\\ & 0,\mathrm{else}.\end{array}$$
The likelihood is $${\left(\frac{1}{\sqrt{2\pi }}\right)}^n(\det \mathrm{\Sigma }{)}^{-\frac{1}{2}}\exp (-\frac{1}{2}(y-m{)}^{\mathrm{T}}{\mathrm{\Sigma }}^{-1}(y-m)),$$ where ${\displaystyle y=(y_1,\cdots ,y_n{)}^{\mathrm{T}}}$. This is a function of ${\displaystyle \mu ,\theta ,\sigma }$, which can be estimated by maximizing the logarithm of the likelihood. ${\displaystyle {\mathrm{\Sigma }}^{-1}}$ makes this computationally expensive. We should use some approximation to skip inverse.
An alternative approach is to find a connection to AR models. We can convert to $${\epsilon }_t=\frac{1}{\theta (B)}(y_t-\mu )=(1-\theta B+{\theta }^2{B}^2-{\theta }^3{B}^3+\cdots )(y_t-\mu ),$$ so that $$y_t-\theta y_{t-1}+{\theta }^2{y}_{t-2}-{\theta }^3{y}_{t-3}+\cdots =\frac{\mu }{1+\theta }+{\epsilon }_t.$$ This requires ${\displaystyle |\theta |<1}$. For this AR model, the likelihood is $${\left(\frac{1}{\sqrt{2\pi }\sigma }\right)}^n\exp (-\frac{1}{2{\sigma }^2}\sum _{t=1}^n{(y_t-\frac{\mu }{1+\theta }-\theta y_{t-1}+{\theta }^2{y}_{t-2}-{\theta }^3{y}_{t-3}+\cdots )}^2).$$ This involves ${\displaystyle y_0,y_{-1},y_{-2},\cdots }$ for which we have no data. We can simply let them to be ${\displaystyle 0}$. Now it becomes $${\left(\frac{1}{\sqrt{2\pi }\sigma }\right)}^n\exp (-\frac{S(\mu ,\theta )}{2{\sigma }^2}),$$ where $$\begin{array}{rl}S(\mu ,\theta )=& {(y_1-\frac{\mu }{1+\theta })}^2+{(y_2-\frac{\mu }{1+\theta }-\theta y_1)}^2+\cdots \\ & +{(y_n-\frac{\mu }{1+\theta }-\theta y_{n-1}+{\theta }^2{y}_{n-2}-\cdots +(-1{)}^{n-1}{\theta }^{n-1}{y}_1)}^2.\end{array}$$
The MLE of ${\displaystyle \mu ,\theta }$ comes from ${\displaystyle \underset{\hat{\mu },\hat{\theta }}{\mathrm{minimize}}S(\mu ,\theta ).}$
This is a nonlinear minimization that can be done in packages in Python like scipy. It's easy to see that ${\displaystyle \hat{\sigma }=\frac{S(\hat{\mu },\hat{\theta })}{n}}$.

For uncertainty quantification, we can take a Bayesian approach. First assume prior $$\theta \sim \mathrm{Uniform}(-1,1),\mu \sim \mathrm{Uniform}(-C,C),\log \sigma \sim \mathrm{Uniform}(-C,C)$$ for a large ${\displaystyle C\to \mathrm{\infty }}$. Note that we restrict ${\displaystyle |\theta |<1}$.
The posterior is then $$\begin{array}{rl}{f}_{\mu ,\theta ,\sigma |\mathrm{data}}(\mu ,\theta ,\sigma )& \propto {\left(\frac{1}{\sqrt{2\pi }\sigma }\right)}^n\exp (-\frac{S(\mu ,\theta )}{2{\sigma }^2})\times \frac{1}{\sigma }\mathbf{1}\{-1<\theta <1,-C<\mu ,\log \sigma <C\}\\ & \propto n^{-n-1}\exp (-\frac{S(\mu ,\theta )}{2{\sigma }^2})\mathbf{1}\{-1<\theta <1,-C<\mu ,\log \sigma <C\}.\end{array}$$ To obtain the posterior of ${\displaystyle \mu ,\theta }$ alone, we integrate the above w.r.t ${\displaystyle \sigma }$. Then we have $$f_{\mu ,\theta |\mathrm{data}}(\mu ,\theta )\propto {\left(\frac{1}{S(\mu ,\theta )}\right)}^{\frac{n}{2}}\mathbf{1}\{-1<\theta <1,-C<\mu <C\}.$$
This can be evaluated numerically over a grid of ${\displaystyle \mu ,\theta }$ and then approximated. Or, we can approximate with a suitable t distribution by doing a Taylor expansion of ${\displaystyle S(\mu ,\theta )}$ near ${\displaystyle \hat{\mu },\hat{\theta }}$. I.e., let ${\displaystyle \alpha =(\mu ,\theta )}$ and ${\displaystyle \hat{\alpha }=(\hat{\mu },\hat{\theta })}$: $$\begin{array}{rl}S(\alpha )& =S(\hat{\alpha })+⟨\mathrm{\nabla }S(\hat{\alpha }),\alpha -\hat{\alpha }⟩+(\alpha -\hat{\alpha }{)}^{\mathrm{T}}(\frac{1}{2}HS(\hat{\alpha }))(\alpha -\hat{\alpha })\\ & =S(\hat{\alpha })+(\alpha -\hat{\alpha }{)}^{\mathrm{T}}(\frac{1}{2}HS(\hat{\alpha }))(\alpha -\hat{\alpha }),\end{array}$$ where we used ${\displaystyle \mathrm{\nabla }S(\hat{\alpha })=0}$ because ${\displaystyle \hat{\alpha }}$ is a minimizer, and ${\displaystyle HS(\hat{\alpha })}$ is the Hessian of ${\displaystyle S}$. Therefore$$\begin{array}{rl}& f_{\mu ,\theta |\mathrm{data}}(\mu ,\theta )\\ \propto & {\left(\frac{1}{S(\mu ,\theta )}\right)}^{\frac{n}{2}}\mathbf{1}\{-1<\theta <1,-C<\mu <C\}\\ \propto & {\left(\frac{S(\hat{\alpha })}{S(\alpha )}\right)}^{\frac{n}{2}}\mathbf{1}\{-1<\theta <1,-C<\mu <C\}\\ =& {\left(\frac{S(\hat{\alpha })}{S(\hat{\alpha })+(\alpha -\hat{\alpha }{)}^{\mathrm{T}}(\frac{1}{2}HS(\hat{\alpha }))(\alpha -\hat{\alpha })}\right)}^{\frac{n}{2}}\mathbf{1}\{-1<\theta <1,-C<\mu <C\}\\ =& {\left(\frac{1}{1+(\alpha -\hat{\alpha }{)}^{\mathrm{T}}(\frac{1}{2S(\hat{\alpha })}HS(\hat{\alpha }))(\alpha -\hat{\alpha })}\right)}^{\frac{n}{2}}\mathbf{1}\{-1<\theta <1,-C<\mu <C\}\\ =& {\left(\frac{1}{1+\frac{1}{n-2}(\alpha -\hat{\alpha }{)}^{\mathrm{T}}(\frac{n-2}{2S(\hat{\alpha })}HS(\hat{\alpha }))(\alpha -\hat{\alpha })}\right)}^{\frac{n-2+2}{2}}\mathbf{1}\{-1<\theta <1,-C<\mu <C\}\end{array}$$
Comparing with $${\left(\frac{1}{1+\frac{1}{k}(x-m{)}^{\mathrm{T}}{\mathrm{\Sigma }}^{-1}(x-m)}\right)}^{\frac{k+p}{2}}$$ for the p-variate t-density ${\displaystyle t_{k,p}(\mu ,\mathrm{\Sigma })}$, we see that $$\alpha |\mathrm{data}\sim t_{n-2,2}(\hat{\alpha },\frac{S(\hat{\alpha })}{n-2}{(\frac{1}{2}HS(\hat{\alpha }))}^{-1}).$$